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ABSTRACT 

A re-analysis of Gliese 581 HARPS and HIRES precision radial velocity data was 
carried out with a Bayesian multi-planet Kepler periodogram (from 1 to 6 planets) 
based on a fusion Markov chain Monte Carlo algorithm. In all cases the analysis 
included an unknown parameterized stellar jitter noise term. For the HARPS data set 
the most probable number of planetary signals detected is 5 with a Bayesian false alarm 
probability of 0.01. These include the 3.1498±0.0005, 5.3687±0.0002, 12.927l^;^J^^, and 
66.9±0.2d periods reported previously plus a 399 j^ijgd period. The orbital eccentricities 

are O.Otj^g, O.OOt^'^Q, QAQt^^, 0.33^;^-?;^, and 0.02l:°-^^, respectively. The semi-major 
axis and Msini of the 5 planets are (0.0285±0.0006 au, 1.9±O.3M0), (0.0406±0.0009 
au, 15.7±O.7M0), (0.073 ±0.002 au, 5.3±O.4M0), (0.218 ±0.005 au, 6.7±0.8Mq), 
and (0.7 ±0.2 au, 6.6t2;?M0), respectively. 

The analysis of the HIRES data set yielded a reliable detection of only the 
strongest 5.37 and 12.9 day periods. The analysis of the combined HIRES/HARPS 
data again only reliably detected the 5.37 and 12. 9d periods. Detection of 4 planetary 
signals with periods of 3.15, 5.37, 12.9, and 66. 9d was only achieved by including an 
additional unknown but parameterized Gaussian error term added in quadrature to 
the HIRES quoted errors. The marginal distribution for the sigma of this additional 
error term has a well defined peak at 1.8 ± 0.4m s~^. It is possible that this additional 
error arises from unidentified systematic effects. We did not find clear evidence for a 
fifth planetary signal in the combined HIRES/HARPS data set. 

Key words: stars: planetary systems; methods: statistical; methods: data analysis; 
techniques: radial velocities. 
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1 INTRODUCTION 

A remarkable array of new ground based and space 
based astronomical tools have finally provided astronomers 
access to other solar systems with over 500 planets 
dis covered to date, starting f rom th e pioneering work 




Campbefi. Walker fc Yand ll 19881) . |Wolszczan & Frail' 
iMavor fc Queloj (jigoi K and Marcy & Butler 
Recent interest has focused on the Gl 581 plane- 
tary system (also designated GJ 581 in the literature). It 
was already known to harbor three planets, including two 
super-Earth planets that straddle its habitable zone: Gl 581 
b with a period of 5.37d (Bonfils et al. 2005b' ) Gl 5 81 c (pe- 
riod 12.9d) and d (period 82 d) (lUdrv et al.jl2007l ). Armed 
with additional HARPS data, iMavor et al.l l|2009l ') reported 
the detection of an additional planet Gl 581e with a mini- 
mum mass of 1.9M0 and a period of 3.15d. They also cor- 
rected previous confusion about the orbital period of Gl 581d 
(the outermost planet) with a one-year alias at 82 days. The 
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revised period is 66. 8d, and positions the semi- major axis 
inside the habitable zone of the low mass star. IVogt et akl 
(|2010| ) reported the analysis of the combined HIRES and 
HARPS data set spanning lly, claiming the detection of two 
additional planets Gl 581f & g. Gl 581f has a period of 433d, 
a minimum-mass of T.OMq, and a semi-major axis of 0.758 
au. Gl 581g has a period of 36. 6d, a minimum-mass 3.1M0, 
and a semi-major axis of 0.146 au. The estimated equilib- 
rium temperature of Gl 581g is 228 K, placing it squarely in 
the middle of the habitable zone of the star. The lVogt et al.l 
(,2010. ) analysis assumed circular orbits for all 6 planets. 



lGregorvll2009l and [Gregory fc Fischeij[2010l presented a 
Bayesian hybrid or fusion MCMC algorithm that incorpo- 
rates parallel tempering (PT), simulated annealing and a 
genetic crossover operation to facilitate the detection of a 
global minimum in . This enables the efficient exploration 
of a large model parameter space starting from a random 
location. When implemented with a multi-planet Kepler 
model, it is able to identify any significant periodic signal 
component in the data that satisfies Kepler's laws and is able 
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to function as a multi-planet Kepler periodogram0. In addi- 
tion, the Bayesian MCMC algorithm provides full marginal 
parameters distributions. The algorithm includes an innova- 
tive adaptive control system that automates the selection of 
efficient parameter propos al distributions even if the param- 
eters are highly correla ted lGregorvll2010l. A rece nt applica- 
tion of the algorithm iGregorv fc Fischer! [2OI0I) confirmed 
the e xistence of a disputed second planet ( Fischer et al.l 
120021 ) in 47 Ursae Majoris (47 UMa) and provided orbital 
constraints on a possible additional long period planet with 
a period ~ lOOOOd. 

This paper reports the results of a re-ana lysis of the 
HARPS (|Mavor et al.l |2009D and HIRES data (|Vogt et all 
I2OIOI ) for Gl 581 using the above Bayesian multi-planet Ke- 
pler periodogram. Section [2] provides an introduction to our 
Bayesian approach and describes the adaptive fusion MCMC 
algorithm. Section [3] gives the model equations and pri- 
ors. Sections 13] and [S] present the parameter estimation and 
model selection results for the analysis of the HARPS data 
alone. Section[6]is devoted to the analysis of the HIRES data 
followed by the analysis of the combination of HIRES and 
HARPS data. The final two sections are devoted to discus- 
sion and conclusions. 



2 THE ADAPTIVE FUSION MCMC 

The adaptive fusion MCMC (FMCMC) is a very general 
Bayesian nonlinear model fitting program. After specifying 
the model, Mi, the data, D, and priors, /, Bayes theo- 
rem dictates the target joint probability distribution for the 
model parameters which is given by 



p{X\D,Nh,I) = C p{X\Kh,I)xp{D\M.,X,I). 



(1) 



where C is the normalization constant and X represent the 
set of model parameters. The first term on the RHS of the 
equation, p{X\Mi, I), is the prior probability distribution of 
X, prior to the consideration of the current data D. The 
second term, p{D\X, Mi, I), is called the likelihood and it is 
the probability that we would have obtained the measured 
data D for this particular choice of parameter vector X, 
model Mi, and prior information I. At the very least, the 
prior information, I, must specify the class of alternative 
models (hypotheses) being considered (hypothesis space of 
interest) and the relationship between the models and the 
data (how to compute the likelihood). In some simple cases 
the log of the likelihood is simply proportional to the familiar 



this type of problem see iGregorvl (|2005bl ). 

To compute the marginals for any subset of the 
parameters it is necessary to integrate the joint probability 
distribution over the remaining parameters. For example, 
the marginal probability density function (PDF) of the 
orbital period in a one planet radial velocity model fit is 



^ Following on from the pion eering work on Bayesian peri- 
odograms bv ljavnej l ll987t) and iBretthorsj lll988l ') 
^ In earlier papers the algorithm was referred to as a hybrid 
MCMC. We subsequently learned that this term already exists in 
the literature in connection with a Hamiltonian version of MCMC. 
In this paper we replace the term hybrid by fusion. 



given by 

p[P\D,Mx,I) = J dK J dV J de J dx J duj J ds 
X p{P,K,V,e,x,uj,s\D,M^,I) 
oc p{P\Mi,I) dK--- I ds 



X p{K,V,e,x,uj,s\Mi,I) 
xpiD\Mi,P,K,V,e,x,co,s,I), (2) 

where p{P, K, V, e, x, 1^, sl^*. Mi, /) is the target joint prob- 
ability distribution of the radial velocity model parameters 
(P, K, V, e, X, 1^) and s is an extra noise parameter which is 
discussed in Section[3l p{P\Mi, I) is the prior for the orbital 
period parameter, p{K,V,e,x,iJ-^,s\Mi,I) is the joint prior 
for the other parameters, and p{D\Mi, P, K,V, e,x,'^, s, I) 
is the likelihood. For a five planet model fit we need to in- 
tegrate over 26 parameters to obtain p{P\D, Mi, I). Inte- 
gration is more difficult than maximization, however, the 
Bayesian solution provides the most accurate information 
about the parameter errors and correlations without the 
need for any additional calculations, i.e., Monte Carlo sim- 
ulations. Bayesian model selection requires integrating over 
all the model parameters. 

In high dimensions, the principle tool for carrying out 
the integrals is Markov chain Monte Carlo based on the 
Metropolis algorithm. The greater efficiency of an MCMC 
stems from its ability, after an initial burn-in period, to gen- 
erate samples in parameter space in direct proportion to the 
joint target probability distribution. In contrast, straight 
Monte Carlo integration randomly samples the parameter 
space and wastes most of its time sampling regions of very 
low probability. 

MCMC algorithms avoid the requirement for com- 
pletely independent samples, by constructing a kind of ran- 
dom walk in the model parameter space such that the num- 
ber of samples in a particular region of this space is pro- 
portional to a target posterior density for that region. The 
random walk is accomplished using a Markov chain, whereby 
the new sample, Xt+i, depends on previous sample Xt ac- 
cording to a time independent entity called the transition 
kernel, p{Xt+i\Xt) ■ The remarkable property of p{Xt+i\Xt) 
is that after an initial burn-in period (which is discarded) 
it generates samples of X with a probability density pro- 
portional to the desired posterior p{X\D, Mi, I) (e.g., see 
Chapter 12 of Gregory 2005 for details). 

The transition kernel, p{Xt+i\Xt) is given by 



p{Xt+i\Xt) = qiXt+i\Xt)a{Xt,Xt+i), 



(3) 



where a{Xt, Xt+i) is called the acceptance probability and 
is given by equation 3] This is achieved by proposing a new 
sample Xt+i from a proposal distribution, q{Xt+i\Xt), which 
is easy to evaluate and is centered on the current sample 
Xt- The proposal distribution can have almost any form. 
A common choice for q{Xt+i\Xt) is a multivariate normal 
(Gaussian) distribution. With such a proposal distribution, 
the probability density decreases with distance away from 
the current sample. The new sample Xt+i is accepted with 
a probability a{Xt, Xt+i) given by 
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a{Xt,Xt+i) = min 



1 p{Xt+i 



D,I) q{Xt\Xt, 



p{Xt\D,I) q{Xt+^\Xt) 



(4) 



where q{Xt\Xt+i) = q{Xt+i\Xt) for a symmetrical proposal 
distribution. If the proposal is not accepted the current sam- 
ple Xt is repeated. 

An important feature that prevents the fusion MCMC 
from becoming stuck in a l ocal p robability maximum 
is parallel tempering (iGeveJ l|l99ll ) and re-invented by 
iHukushima fc Nemotd (| 19961 ) under the name exchange 
Monte Carlo). Multiple MCMC chains are run in parallel. 
The joint distribution for the parameters (X) of model Mi, 
for a particular chain, is given by 



Tv{X\D,M,,I,l3) cxp{X\M,,I) xp{D\X,M,,lf 



(5) 



Each MCMC chain corresponding to a different jS, with the 
value of P ranging from zero to 1. When the exponent /? = 1, 
the term on the LHS of the equation is the target joint prob- 
ability distribution for the model parameters, p{X\D, Mi, I). 
For -C 1, the distribution is much flatter. 

In equation [S] an exponent l3 — yields a joint den- 
sity distribution equal to the prior. The reciprocal of /? is 
analogous to a temperature, the higher the temperature the 
broader the distribution. For parameter estimation purposes 
8 chains with 

P = {0.09,0.13,0.20,0.29,0.39,0.52,0.72,1.0} were em- 
ployed. At an interval of 10 iterations, a pair of adjacent 
chains on the tempering ladder are chosen at random and 
a proposal made to swap their parameter states. A Monte 
Carlo acceptance rule determines the probability for the pro- 
posed swap to occur (e.g., Gregory 2005a, equation 12.12). 
This swap allows for an exchange of information across the 
population of parallel simulations. In low /3 (higher tem- 
perature) simulations, radically different configurations can 
arise, whereas in higher /? (lower temperature) states, a con- 
figuration is given the chance to refine itself. The lower /3 
chains can be likened to a series of scouts that explore the 
parameter terrain on different scales. The final samples are 
drawn from the P — 1 chain, which corresponds to the de- 
sired target probability distribution. The choice of /3 val- 
ues can be checked by computing the swap acceptance rate. 
When they are too far apart the swap rate drops to very low 
values. A swap acceptance rate of « 40% works well. 

At each iteration, a single joint proposal to jump to 
a new location in the parameter space is generated from 
independent Gaussian proposal distributions (centered on 
the current parameter location), one for each parameter. In 
general, the cr's of these Gaussian proposal distributions are 
different because the parameters can be very different enti- 
ties. If the (t's are chosen too small, successive samples will 
be highly correlated and will require many iterations to ob- 
tain an equilibrium set of samples. If the cr's are too large, 
then proposed samples will very rarely be accepted. The 
process of choosing a set of useful proposal cr's when deal- 
ing with a large number of different parameters can be very 
time consuming. In parallel tempering MCMC, this prob- 
lem is compounded because of the need for a separate set of 
Gaussian proposal a's for each tempering chain. This pro- 
cess is automat ed by an innovative three sta ge statistical 
control system l|Gregorvll2007bl . lGregorvl[2009l ) in which the 
error signal is proportional to the difference between the 
current joint parameter acceptance rate and a target accep- 



tance rate, typically 25% (|Roberts et al.|[l997t ). A schematic 
of the first two stages of the adaptive control system (CS) is 
shown in Fig. 1. A third stage that handles highly correlated 
parameters is described in Section [2. II 

The first stage CS, which involves annealing the set of 
Gaussian proposal distribution a's, was described in Gre- 
gory 2005a. An initial set of proposal cr's (~ 10% of the 
prior range for each parameter) are used for each chain. Dur- 
ing the major cycles, the joint acceptance rate is measured 
based on the current proposal cr's and compared to a target 
acceptance rate. During the minor cycles, each proposal a 
is separately perturbed to determine an approximate gradi- 
ent in the acceptance rate. The cr's are then jointly modified 
by a small increment in the direction of this gradient. This 
is done for each of the parallel chains. Proposals to swap 
parameter values between chains are allowed during major 
cycles but not within minor cycles. 

The annealing of the proposal cr's occurs while the FM- 
CMC is homing in on any significant peaks in the target 
probability distribution. Concurrent with this, another as- 
pect of the annealing operation takes place whenever the 
Markov chain is started from a location in parameter space 
that is far from the best fit values. This automatically arises 
because all the m odels considered incorporate an extra ad- 
ditive noise term (|Gregorvll2005bl ) whose probability distri- 
bution is Gaussian with zero mean and with an unknown 
standard deviation s. When the of the fit is very large, 
the Bayesian Markov chain automatically inflates s to in- 
clude anything in the data that cannot be accounted for by 
the model with the current set of parameters and the known 
measurement errors. This results in a smoothing out of the 
detai l ed str ucture in the surface and, as pointed out by 
iFordI I 20061 ). allows the Markov chain to explore the large 
scale structure in parameter space more quickly. The chain 
begins to decrease the value of the extra noise as it settles in 
near the best-fit parameters. An example of this is shown in 
Fig. [5] for a two planet fit to the HARPS data as discussed 
in Section [4.11 The three panels shows the evolution of the 
Logio [Prior x Likelihood] , the s parameter and the two pe- 
riod parameters. In the early stages the extra noise is inflated 
to around 16m s~^ and then rapidly decays to much lower 
values as it homes in on possible solutions. This is similar to 
simulated annealing, but does not require choosing a cooling 
scheme. In this example, the starting parameter values were 
far from the best and the MCMC algorithm finds several less 
probable solutions on route to a final best choice. Initially it 
homes in on the two periods of 5.37 and 66. 9d and the CS 
switches off around iteration 220,000, but around iteration 
500,000 it detects a much more probable solution for the 
period combination of 5.37 and 12. 9d. In this example the 
adaptive control system switched on again briefly following 
the detection of the much improved solution. 

Although the first stage CS achieves the desired joint 
acceptance rate, it often happens that a subset of the pro- 
posal cr's are too small leading to an excessive autocorrela- 
tion in the FMCMC iterations for these parameters. Part of 
the second stage CS corrects for this. The goal of the sec- 
ond stage is to achieve a set of proposal cr's that equalizes 
the FMCMC acceptance rates when new parameter values 
are proposed separately and achieves the desired acceptance 
rate when they are pr oposed jointly. D etails of the second 
stage CS were given in lGregorvll2007bl . 
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Figure 1. The first of two scliematics on tlie operation of the adaptive fusion MCMC (FMCMC) algorithm. 



The first stage is run only once at the beginning, but the 
second stage can be executed repeatedly, whenever a signifi- 
cantly improved parameter solution emerges. Frequently, the 
algorithm homes in on the most significant peak within the 
span of the first stage CS and the second stage improves the 
choice of proposal cr's based on the highest probability pa- 
rameter set. Occasionally, a new higher (by a user specified 
threshold) target probability parameter set emerges after 
the first two stages of the CS have completed. The control 
system has the ability to detect this and automatically re- 
activate the second stage. In this sense the CS is adaptive. 
If this happens the iteration corresponding to the end of 
the control system is reset. The requirement that the transi- 
tion kernel be time independent means that q{Xt+i\Xt) be 
time independent, so useful FMCMC simulation data are 
obtained only after the CS is switched off. 

The adaptive capability of the control system can be 
appreciated from an examination of Fig. 1. The upper left 
portion of the figure depicts the FMCMC iterations from the 
8 parallel chains, each corresponding to a different tempering 
level P as indicated on the extreme left. One of the outputs 
obtained from each chain at every iteration (shown at the 
far right) is the log prior + log likelihood. This information 
is continuously fed to the CS which constantly updates the 
most probable parameter combination regardless of which 
chain the parameter set occurred in. This is passed to the 
'Peak parameter set' block of the CS. Its job is to decide if a 
significantly more probable parameter set has emerged since 
the last execution of the second stage CS. If so, the second 
stage CS is re-run using the new more probable parameter 
set which is the basic adaptive feature of the existing CS. 



The CS also includes a genetic algorithm block which is 
shown in the bottom right of Fig. 1. The current parameter 
set can be treated as a set of genes. In the present version, 
one gene consists of the parameter set that specifies one or- 
bit. On this basis, a three planet model has three genes. At 
any iteration there exist within the CS the most probable 
parameter set to date Xmax, and the current most proba- 
ble parameter set of the 8 chains, Xcui- At regular intervals 
(user specified) each gene from Xcm is swapped for the cor- 
responding gene in Xmax- If either substitution leads to a 
higher probability it is retained and Xmax updated. The ef- 
fectiveness of this operation can be tested by comparing the 
number of times the gene crossover operation gives rise to 
a new value of Xmax compared to the number of new Xmax 
arising from the normal parallel tempering FMCMC itera- 
tions. The gene crossover operations prove to be very effec- 
tive, and give rise to new Xmax values ~ 3 times more often. 
Of course, most of these swaps lead to very minor changes 
in probability but occasionally big jumps are created. 



Gene swaps from Xcur2, the parameters of the second 
most probable current chain, to Xmax can also be also uti- 
lized. This gives rise to new values of ^max at a rate approx- 
imately half that of swaps from Xcur to Xmax. Crossover op- 
erations at a random point in the entire parameter set did 
not prove as effective except in the single planet case where 
there is only one gene. Further experimentation with this 
concept is ongoing. 
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Figure 3. Schematic outlining the operation of the third stage of the adaptive fusion MCMC algorithm that handles correlated param- 
eters. 



2.1 Highly correlated parameters 

The part of the algorithm described above (Figure 1) is most 
efhcient when working with model parameters that are in- 
dependent or transformed to new independent parameters. 
New parameter values are jointly proposed based on inde- 
pendent Gaussian proposal distributions (T scheme), one 
for each parameter. Initially, only this T proposal system is 
used and it is clear that if there are strong correlations be- 
tween any parameters the cr's of the independent Gaussian 
proposals will need to be very small for any proposal to be 
accepted and consequently convergence will be very slow. 
However, the accepted T proposals will generally cluster 
along the correlation path. In the optional third stage of the 
control system shown in Figure im every second accepted T 
proposal is appended to a correlated sample buffer. Only the 
300 most recent additions to the buffer are retained. A 'C' 
proposal is generated from the difference between a pair of 
randomly selected samples drawn from the correlated sam- 
ple buffer, after multiplication by a constant. The value of 
this constant is computed automatically by another control 
system module which ensures that the 'C' proposal accep- 
tance rate is close to 25%. With very little computational 
overhead, the 'C' proposals provide the scale and direction 
for efficient jumps in a correlated parameter space. 

The final proposal distribution is a random selection of 
T and 'C' proposals such that each is employed 50% of the 
time. The overhead to generate the 'C' proposals is minimal. 
The combination ensures that the whole parameter space 
can be reached and that the FMCMC chain is aperiodic. 



The parallel tempering feature operates as before to avoid 
becoming trapped in a local probability maximum. 

Because the 'C' proposals reflect the parameter correla- 
tions, large jumps are possible allowing for much more effl- 
cient movement in parameter space than can be achieved by 
the T' proposals alone. Once the first two stages of the con- 
trol system have been turned off, the third stage continues 
until a minimum of an additional 300 accepted T proposals 
have been added to the buffer and the 'C' proposal accep- 
tance rate is within the range ^ 0.22 and ^ 0.28. At this 
point further additions to the buffer are terminated and this 
sets a lower bound on the burn-in period. 

Full details on the operation and testi ng of the com- 
bined T and 'C' proposal scheme are given in lGregorvll2010l . 



3 MODELS AND PRIORS 

In this section we describe the model fitting equations and 
the selection of priors for the model parameters. We have 
investigated the Gl 581 data using models ranging from 1 
to 6 planets. For a one planet model the predicted radial 
velocity is given by 

viU) = V + K[cos{9(ti + xP) + uj} + ecosuj], (6) 

and involves the 6 unknown parameters 

V = a, constant velocity. 

K — velocity semi-amplitude. 

P = the orbital period. 

e = the orbital eccentricity. 
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Figure 2. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus MCMC iteration for a blind 2 planet fit to the 
HARPS data. The middle panel is a similar plot for the extra 
noise term s. Initially s is inflated and then rapidly decays to 
a much lower level as the better fit parameter values are ap- 
proached. The lower panel shows the values of the two unknown 
period parameters versus iteration number. The two starting pe- 
riods of 2.5 and 20d are shown on the left hand side of the plot 
at a negative iteration number. 



u) = the longitude of periastron. 

X = the fraction of an orbit, prior to the start of data 
taking, that periastron occurred at. Thus, xP — the number 
of days prior to ti = that the star was at periastron, for 
an orbital period of P days. 

6{ti + xP) = the true anomaly, the angle of the star in 
its orbit relative to periastron at time ti. 

We utilize this form of the equation because we obtain 
the dependence of 6 on ti by solving the conservation of 
angular momentum equation 

de_2n[l + e cos e{t^ + x P)? 
dt 



P(l 



- 0. 



(7) 



Our algorithm is implemented in Mathematica and it proves 
faster for Mathematica to solve this differential equation 
than solve the equations relating the true anomaly to the 
mean anomaly via the eccentric anomaly. Mathematica gen- 
erates an accurate interpolating function between t and 9 
so the differential equation does not need to be solved sepa- 



rately for each ti. Evaluating the interpolating function for 
each ti is very fast compared to solving the differential equa- 
tion, so the algorithm should be able to handle much larger 
samples of radial velocity data than those currently available 
without a significant increase in computational time. For ex- 
ample, an increase in the data by a factor of 6.5 resulted in 

only an 18% increase in execution t ime. 

As described in more detail in lGregorvl[2007al , we em- 
ployed a re-parameterization of x ^-nd uj to improve the 
MCMC convergence speed motivated by the work of Ford 
(2006). The two new parameters are ^ = 27rx + ^ and 
(f> = 2-KX — ^. Parameter is well determined for all ec- 
centricities. Although (j) is not well determined for low ec- 
centricities, it is at least orthogonal to the t/; parameter. We 
use a uniform prior for i/; in the interval to 47r and uniform 
prior for (ft in the interval — 27r to -|-27r. This insures that a 
prior that is wraparound continuous in (x, i^) maps into a 
wraparound continuous distribution in ("i/),!^). To account 
for the Jacobian of this re-parameterization it is necessary 
to multiply the Bayesian integrals by a factor of (47r) 
where nplan = the number of planets in the model. Also, 
by utilizing the orthogonal combination [il^,(j>) it was not 
necessary to make use of the 'C proposal scheme outlined 
in Section [2.11 which typically saves about 25% in execution 
time. 

In a Bayesian analysis we need to specify a suitable prior 
for each parameter. These are tabulated in Table [T] For the 
current problem, the prior given in Equation[5]is the product 
of the individual parameter priors. Detailed argume nts for 
the choice of each p rior were given in lGregorvl[2007al . 

iGregorvl l2007al discussed two different strategies to 
search the orbital frequency parameter space for a multi- 
planet model: (i) an upper bound on /i ^ /2 ^ • ■ • ^ /n 
is utilized to maintain the identity of the frequencies, and 
(ii) all fi are allowed to roam over the entire frequency 
range and the parameters re-labeled afterwards. Case (ii) 
was found to be significantly more successful at converging 
on the highest posterior probability peak in fewer iterations 
during repeated blind frequency searches. In addition, case 
(ii) more easily permits the identification of two planets in 
1:1 resonant orbits. We adopted approach (ii) in the current 
analysis. 

All of the models considered in this paper incorporate 
an extra noise parameter, s, that can allow for any addi- 
tional noise beyond the known measurement uncertainties0. 
We assume the noise variance is finite and adopt a Gaussian 
distribution with a variance s^. Thus, the combination of 
the known errors and extra noise has a Gaussian distribution 
with variance = a'l+s^ , where Oi is the standard deviation of 
the known noise for i"* data point. For example, suppose that 
the star actually has two planets, and the model assumes 
only one is present. In regard to the single planet model, the 
velocity variations induced by the unknown second planet 
acts like an additional unknown noise term. Other factors 
like star spots and chromospheric activity can also con- 



^ In the absence of detailed knowledge of the sampling distribu- 
tion for the extra noise, we pick a Gaussian because for any given 
finite noise variance it is the distribution with the largest uncer- 
tainty as me asured by the entropy, i.e., the maximum entropy 
distribution ll Javneslll957l , lGregorvll2005al section 8.7.4.) 
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Table 1. Prior parameter probability distributions. 



Parameter 



Prior 



Lower bound Upper bound 



Orbital frequency 

Velocity Ki 
(m s-i) 

V (m s-l) 

Si Eccentricity 

X orbit fraction 

{jji Longitude of 
periastron 



p(ln /i , In /2, • ■ • In /„ | Af„, /) = l^n^f^jf^^^. 
(n =number of planets) 

Modified Jeffreys^ 



Uniform 

a) Uniform 

b) Ecc. noise bias correction filter 
Uniform 

Uniform 



Extra noise (m s -*-) — (^+^0) 



1/1.1 d 1/1000 yr 



(Ko = l) Xmax {^Y^^ 



K^B.^ = 2129 



1 

0.99 
1 

27r 



In 1 + ; 



(so = 1) Kn 



°' Since the prior lower limits for K and s include zero, we used a modified Jeffreys prior of the form 

p{X\M,I) = i 2— (8) 

X + Xo ln(l + 2^) 

For X <C Xq, p{X\M, I) behaves like a uniform prior and for X ^ Xo it behaves like a Jelfreys prior. The 
In (1 + '^Xq'^ ) term in the denominator ensures that the prior is normalized in the interval to Xmax. 



tribute to this extra velocity noise term which is often re- 
ferred to as stellar jitter. Several researchers have attempted 
to estimate stellar jitter for individual st ars based on sta- 
tistical correlations wi t h observables (e.g.. [Saar fc Donahue! 
I1997I . ISaar et al.iri998l . I Wrightl [2OO5I ) . In general, nature is 
more complicated than our model and known noise terms. 
Marginalizing s has the desirable effect of treating anything 
in the data that can't be explained by the model and known 
measurement errors as noise, leading to conservative esti- 
mates of orbital parameters. See Sections 9.2.3 and 9.2.4 of 
iGregorvl (|2005al ) for a tutorial demonstration of this point. 
If there is no extra noise then the posterior probability dis- 
tribution for s will peak at s = 0. The upper limit on s was 
set equal to i^max- We employed a modified Jeffrey's prior 
for s with a knee, so = Im s~^. 

We used two different choices of priors for eccentricity, 
a uniform prior and eccentricity noise bias correction filter 
that is described in the next section. 



3.1 Eccentricity bias 

In lGreeorv fc Fischer|[2O10l . the velocities of model fit resid- 
uals were randomized in multiple trials and processed using 
a one planet version of the FMCMC Kepler periodogram. 
In this situation periodogram probability peaks are purely 
the result of the effective noise. The orbits corresponding 
to these noise induced periodogram peaks exhibited a well 



defined statistical bias towards high eccentricity. They of- 
fered the following explanation for this effect. To mimic a 
circular velocity orbit the noise points need to be correlated 
over a larger fraction of the orbit than they do to mimic a 
highly eccentric orbit. For this reason it is more likely that 
noise will give rise to spurious highly eccentric orbits than 

low eccentricity orbits. 

[Gregory fc Fischerll2O10l characterized this eccentricity 
bias and designed a correction filter that can be used as 
an alternate prior for eccentricity to enhance the detection 
of planetary orbits of low or moderate eccentricity. On the 
basis of our understanding of the mechanism underlying the 
eccentricity bias, we expect the eccentricity prior filter to 
be generally applicable to searches for low amplitude orbital 
signals in precision radial velocity data sets. The probability 
density function for this filter is shown by the solid black 
curve in Figure |4] and is given by 

pdf{e) = 1.3889-1.5212eV0.53944e^-1.6605(e-0.24821)^(9) 

In a related study, IShen and Turner! (|2009l ) explored least- 
Keplerian fits to radial velocity data using synthetic data 
sets. They found that the best fit eccentricities for low signal- 
to-noise ratio K/a ^ 3 and moderate number of observa- 
tions Nobs ^ 60, were systematically biased to higher values, 
leading to a suppression of the number of nearly circular or- 
bits. 

In the analysis of the Gl 581 data we used the eccen- 
tricity noise bias correction filter as the eccentricity prior on 
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0.4 0.6 

Eccentricity 



Figure 4. The best fit polynomial (thin black curve) to the 
reciprocal of the mean of the eccentricity bias determined by 
iGregorv fc Fischer! |2010| . After normalization this yields the ec- 
centricity noise bias correction filter (lower solid black curve). 



fits of all the models, with occasional runs using a uniform 
eccentricity prior to test the robustness of our conclusions. 



4 ANALYSIS OF THE HARPS DATA 

The HARPS data (|Mavor et all |2009| ) were retrieved elec- 
tronically 0. A mean velocity of -9.2080205km s~^ was sub- 
tracted and the remainder converted to units of m s~^. 
Panel (a) of Figure [S] hows the HARPS observations of Gl 
581. Panel (b) shows a blow-up of a portion of the mean 5 
planet model fit compared to the data, and panel (c) shows 
the residuals. The zero reference time is the mean time of 
the HARPS observations which corresponds to a Julian day 
number = 2, 454, 186.6178. 
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Figure 5. Panel (a) shows the HARPS observations of Gl 581. 
Panel (b) shows a blow-up of the mean 5 planet model fit com- 
pared to the data, and panel (c) shows the residuals. 



4.1 Two planet model 

The results of our 2 planet Kepler periodogram analysis of 
this data are shown in Figures [2] and [6l The upper panel 
of Figure [2] shows a plot of the Logio [Prior x Likelihood] 
versus FMCMC iteration for a 2 planet fit of the HARPS 
data. The lower panel shows the values of the two unknown 
period parameters versus iteration number. The two starting 
periods of 2.5 and 20d are shown on the left hand side of the 
plot at a negative iteration number. The larger of the two 
period parameter finds both the 67 and 12. 9d periods but 
the latter has a much larger Logio [Prior x Likelihood] value. 
The median value of the extra noise parameter s — 2.37m 
s-i. 

Figure |6] shows plot of a sample of the FMCMC two pe- 
riod parameters versus a normalized value of Logio [Prior x 
Likelihood], i.e., a 2 planet periodogram. Only values within 
18 decades of the maximum Logio [Prior x Likelihood] are 
plotted but without regard to whether the values occurred 
before or after burn-in. The two prominent periods are 5.37 
& 12. 9d. The second period parameter exhibited many other 
peaks but these were all at least 8 decades less probable. 



|http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A-|-A/507/48"7] 



Q. 




1000 3000 10000 



Figure 6. A plot of the two period parameter values versus a 
normalized value of Logio [Prior X Likelihood] for the 2 planet 
FMCMC Koplor periodogram of the HARPS data. 



4.2 Three planet model 

The results of our 3 planet Kepler periodogram analysis are 
shown in Figures [71 El [9l & \W\ The upper panel of Figure [71 
shows a plot of the Logio [Prior x Likelihood] versus FM- 
CMC iteration for a 3 planet fit of the HARPS data. The 
lower panel shows the FMCMC values of the three unknown 
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Figure 7. The upper panel is a plot of the Logio [Prior X Likeli- 
hood] versus iteration for the three planet FMCMC Kepler peri- 
odogram of the HARPS data. The lower shows the values of the 
three unknown period parameters versus iteration number. The 
three starting periods of 2.5, 20, lOOd are shown on the left hand 
side of the plot at a negative iteration number. 



period parameters versus iteration number. The three start- 
ing periods of 2.5, 20, lOOd are shown on the left hand side of 
the plot at a negative iteration number. The burn-in period 
for this run was 0.13 x 10^ iterations. 

Figure |8] shows plot of a sample of the FMCMC 
three period parameters versus a normalized version of 
Logioprior x Likelihood], i.e., a 3 planet periodogram. Only 
values within 5 decades of the maximum Logio [Prior x Like- 
lihood] are plotted but without regard to whether the values 
occurred before or after burn-in. Three prominent periods 
were clearly detected: 5.37, 12.9, 66. 9d. The third period pa- 
rameter exhibited four other peaks but these were all more 
than 2 decades less probable. The most probable of these 
has a period ~ 413d. The spectral peak at 82d coincides 
with a one year alias (1/67— 1/365 ~ 1 /82) of the dominant 
67d p eriod. For more on RV aliases see lDawson fc Fabrvckvl 
l|2010h . The 59d peak is close but not coincident with the 
other one year alias (1/67-1- 1/365 ~ 1/57). 

Figure |9] shows a plot of eccentricity versus period for a 
sample of the FMCMC parameter samples for the 3 planet 
model. There is clearly a large uncertainty in the eccentricity 
of the 67d period which extends down to low eccentricities. 
The 413d period peak exh ibits very large eccentricity values. 
iGreeorv fc Fischeij ()2010[ ) showed that it is more likely that 
noise will give rise to spurious highly eccentric orbits than 
low eccentricity orbits. To mimic a circular velocity orbit 
the noise points need to be correlated over a larger fraction 
of the orbit than they do to mimic a highly eccentric or- 
bit. Even though we are using the noise ind uced eccentricity 
prior proposed in lCreeorv fc Fischeij ()2010l ) we still observe 
a preponderance of high eccentricity orbital solutions in the 
low K value regime. 



30 5967 82 100 

Period (d) 



Figure 8. A plot of the 3 period parameter values versus a nor- 
malized value of Logio [Prior X Likelihood] for the 3 planet FM- 
CMC Kepler periodogram of the HARPS data. 
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Period (d) 

Figure 9. A plot of eccentricity versus period for the 3 planet 
FMCMC Kepler periodogram of the HARPS data. 



Figure [TO] shows a plot of a subset of the FMCMC pa- 
rameter marginal distributions for the 3 planet fit of the 
HARPS data after filtering out the post burn-in FMCMC 
iterations that correspond to the 3 dominant period peaks at 
5.37, 12.9, and 66. 9d. The bottom panel shows the marginal 
for the unknown standard deviation, s, of the additive Gaus- 
sian extra noise term which has a median value of 1.93 m 
s~^. The Bayesian analysis automatically infiates s to ac- 
count for anything in the data that the model and quoted 
measurement errors cannot account for including stellar jit- 
ter. 

The three planet model was also run using a fiat uni- 
form eccentricity to compare with the results obtained with 
the noise induced eccentricity prior. Figure [TT1 shows a com- 
parison of the eccentricity marginals for the noise induced 
prior (solid black curve) and the uniform prior (grey dashed 



4.3 Four planet model 

A one planet fit to the residuals of the three planet fit above 
yielded a dominant Keplerian orbit with a period of 3.15d. 
The results of our 4 planet Kepler periodogram analysis are 
shown in Figures [12] and 1131 All 4 period parameters were 
free to roam within a search range extending from l.ld to 
10 X the data duration. Another run that extended the pe- 
riod search range down to 0.5d yielded the same 4 periods. 
The median value of extra noise parameter s = 1.36m s~^. 
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Figure 10. A plot of a subset of the FMCMC parameter marginal 
distributions for a 3 planet fit of the HARPS data. 





Figure 11. The eccentricity marginals for (a) the noise induced 
eccentricity prior (solid black curve) and (b) the uniform prior 
(grey dashed curve). 
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Figure 12. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus iteration for the four planet FMCMC Kepler peri- 
odogram of the HARPS data. The lower shows the values of the 
four unknown period parameters versus iteration number. 




Figure 13. A plot of a subset of the FMCMC parameter marginal 
distributions for a 4 planet fit of the HARPS data. 



4.4 Five planet model 

The results of a 5 planet Kepler periodogram analysis of 
the HARPS data are shown in Figures [H ^ and [TH] The 
starting period for the fifth period was set = 300d and the 
most probable period found to be ~ 400d. The best set 
of parameters from the 4 planet fit were used as start pa- 
rameters. The fifth period parameter shows 3 peaks but the 
400d period is almost 1000 times stronger than the other 
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Figure 14. A plot of the 5 period parameter values versus a 
normalized value of Logio [Prior X Likelihood] for the 5 planet 
FMCMC Kepler periodogram of the HARPS data. The fifth pe- 
riod parameter points are shown in orange. 
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Figure 15. A plot of eccentricity versus period for the 5 planet 
FMCMC Kepler periodogram of the HARPS data. The fifth pe- 
riod parameter points are shown in grey. 



two. We filtered the 5 planet MCMC results to only include 
fifth period val u es fro m the 400d peak region and used the 
iGelman-Rubiiil (|l992l ) statistic to test for convergence. In 
parallel tempering MCMC, new widely separated parame- 
ter values are passed up the line to the /3 = 1 simulation and 
are occasionally accepted. Roughly every 100 iterations the 
/3 = 1 simulation accepts a swap proposal from its neigh- 
boring simulation. The final /? = 1 simulation is thus an 
average of a very large number of independent P — 1 sim- 
ulations. We divide the /3 = 1 iterations into 10 equal time 
intervals and inter-compared the 10 different essentially in- 
dependent average distributions for each parameter using a 
Gelman- Rubin test. For the five planet model results the 
Gelman-Rubin statistic was ^ 1.01. 

Figure \W\ shows a plot of a subset of the FMCMC pa- 
rameter marginal distributions for the 5 planet fit of the 
HARPS data after fihering out the post burn-in FMCMC 
iterations that correspond to the 5 dominant period peaks 
at 3.15, 5.37, 12.9, 66.9, and 400d. The median value of the 
extra noise parameter s — 1.16m s~^. 




Figure 16. A plot of a subset of the FMCMC parameter marginal 
distributions for a 5 planet fit of the HARPS data. 



4.5 Six planet model 

We also carried out a 6 planet Kepler periodogram analysis 
of the HARPS data and the results are shown in Figures [TT] 
I18II19I and 1201 The best set of parameters from the 5 planet 
fit were used as start parameters and the starting period 
for sixth period was set = 36d. The most probable sixth 
period found was 34.4d. A 34. 4d period also appeared as a 
secondary peak in the 5 planet fit and is evident in Figure [T4l 
Figure \W\ shows a plot of eccentricity versus period for the 
5 planet FMCMC Kepler periodogram. 

Figure [20] shows a plot of a subset of the FMCMC pa- 
rameter marginal distributions for the 6 planet fit of the 
HARPS data after filtering out the post burn-in FMCMC 
iterations that correspond to the 6 dominant period peaks 
at 3.15, 5.37, 12.9, 34.4, 66.9, and 400d. The median value 
of extra noise parameter s = 1.00m s~^. There is consid- 
erable agreement with the 4 and 5 planet marginals shown 
earlier which leads us to conclude that the 66. 9d orbit is 
significantly eccentric with an e « 0.34. There is an indica- 
tion that e « 0.12 for the 12.9d orbit. If the 34.4d orbit is 
real it would also appear to have a significant eccentricity of 

Phase plots for the 6 planet model are shown in Fig- 
ure 1211 The top left panel shows the data and model fit 
versus 3.15d orbital phase after removing the effects of the 
five other orbital periods. To construct this phase plot we 
first filter out the post burn-in FMCMC iterations that cor- 
respond to the 6 dominant period peaks at 3.15, 5.37, 12.9, 
34.4, 66.9, and 400d. The FMCMC output for each of these 
iterations is a vector of the 6 planet orbital parameter set 
plus V . To compute the 3.15d phase plot data we subtract 
the mean velocity curve for the other 5 planets plus V from 
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Figure 17. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus iteration for the FMCMC six planet fit of the 
HARPS data. The lower shows the values of the six unknown 
period parameters versus iteration number. The six starting pe- 
riods are shown on the left hand side of the plot at a negative 
iteration number. 
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Figure 18. A plot of the 6 period parameter values versus a 
normalized value of Logio [Prior X Likelihood] for the 6 planet 
FMCMC Kepler periodogram of the HARPS data. The fifth pe- 
riod parameter points are shown in orange and the sixth in black. 



the measured set of velocities. This is done by taking a sam- 
ple of typically 200 FMCMC iterations and for each iteration 
we compute the predicted velocity points for that realization 
of the 5 planet plus V parameter set. We then construct the 
average of these model prediction data sets and subtract 
that from the data points. These residuals for the set of ob- 
servation times are converted to residuals versus phase using 
the mode of the marginal distribution for the 3.15d period 
parameter. An orbital phase model velocity fit is then com- 
puted at 100 phase points for each realization of the 3.15d 
planet parameter set obtained in the same sample of 200 
iterations as above. At each of these 100 phase points we 
construct the mean model velocity fit and mean ±1 stan- 
dard deviation. The red and blue solid curves in Figure [211 
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Figure 19. A plot of eccentricity versus period for the 6 planet 
FMCMC Kepler periodogram of the HARPS data. 




Figure 20. A plot of a subset of the FMCMC parameter marginal 
distributions for a 6 planet fit to the HARPS data. 



are the mean FMCMC model fit ±1 standard deviation. 
Thus, 68.3% of the FMCMC model fits fall between these 
two curves. 

The other panels correspond to phase plot for the other 
five periods. In each panel the quoted period is the mode of 
the marginal distribution. It is clear that for the 3.15, 5.37, 
12.9, 66. 9d periods the separation of the fit curves are small 
compared to the amplitude. For the 395d period phase plot, 
the wide range of possible orbits that can fit between the 
red and blue curves is refiected by the broad extent of the 
marginal distributions of the parameters. 
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Figure 21. Phase plots for the 6 planet model fit to the IfARPS 
data. The top left panel shows the data and model fit versus 
3.15 day orbital phase after removing the effects of the five other 
orbital periods. The red and blue curves are the mean FMCMC 
model fit ±1 standard deviation. The other five panels correspond 
to phase plot for the other five periods. 



5 MODEL SELECTION FOR HARPS 
ANALYSIS 

One of the great strengths of Bayesian analysis is the built- 
in Occam's razor. More complicated models contain larger 
numbers of parameters and thus incur a larger Occam 
penalty, which is automatically incorporated in a Bayesian 
model se lection analysis in a quantitative fashion (see for 
example. lGregorvll2005al . p. 45). The analysis yields the rel- 
ative probabihty of each of the models explored. 

To compare the posterior probability of the i**" planet 
model to the four planet model we need to evaluate the 
odds ratio, On = p[Mi\D,I)/p(M4,\D,I), the ratio of the 
posterior probability of model Mi to model M4. Application 
of Bayes theorem leads to. 



On = 



p{M,\I) p{D\M,,I) ^ p{Mi\I) 
p{M4I) p{D\Mi,I) - p[Mi\I) 



(10) 



where the first factor is the prior odds ratio, and the second 
factor is called the Bayes factor, Bu. The Bayes factor is 
the ratio of the marginal (global) likelihoods of the models. 
The marginal likelihood for model Mi is given by 



p{D\Mi,I) 



dXp{X\M„I) xp{D\X,M„I). 



(11) 



Thus Bayesian model selection relies on the ratio of marginal 
likelihoods, not maximum likelihoods. The marginal likeli- 
hood is the weighted average of the conditional likelihood, 
weighted by the prior probability distribution of the model 
parameters and s. This procedure is referred to as marginal- 
ization. 

The marginal likelihood can be expressed as the prod- 
uct of the maxi mum likelihood and the Occam penalty (see 
iGregorvl [2005al . page 48). The Bayes factor will favor the 
more complicated model only if the max;imum likelihood ra- 
tio is large enough to overcome this penalty. In the simple 
case of a single parameter with a uniform prior of width 
AX, and a centrally peaked likelihood function with char- 



acteristic width SX, the Occam factor is « 5X/AX. If the 
data is useful then generally SX <^ AX. For a model with 
m parameters, each parameter will contribute a term to the 
overall Occam penalty. The Occam penalty depends not only 
on the number of parameters but also on the prior range of 
each parameter (prior to the current data set, D), as sym- 
bolized in this simplified discussion by AX. If two models 
have some parameters in common then the prior ranges for 
these parameters will cancel in the calculation of the Bayes 
factor. To make good use of Bayesian model selection, we 
need to fully specify priors that are independent of the cur- 
rent data D. The sensitivity of the marginal likelihood to the 
prior range depends on the shape of the prior and is much 
greater for a u niform prior than a Jeffreys prior (e.g., see 
lGregorvll2005al . page 61). In most instances we are not par- 
ticularly interested in the Occam factor itself, but only in the 
relative probabilities of the competing models as expressed 
by the Bayes factors. Because the Occam factor arises au- 
tomatically in the marginalization procedure, its effect will 
be present in any model selection calculation. Note: no Oc- 
cam factors arise in parameter estimation problems. Param- 
eter estimation can be viewed as model selection where the 
competing models have the same complexity so the Occam 
penalties are identical and cancel out. 

The MCMC algorithm produces samples which are in 
proportion to the posterior probability distribution which 
is fine for parameter estimation but one needs the pro- 
portionalit y constant f o r est imating the model marginal 
hkelihood. IClvde et al.l (|2006l ) reviewed the state of tech- 
niques for model select ion from a statistics perspective and 
iFord fc Gregorvl (|2006l ) have evaluated the performance of 
a variety of marginal likelihood estimators in the exoplanet 
context. 

Estimating the marginal likelihood is a very big chal- 
lenge for models with large numbers of parameters, e.g., our 
six planet model has 32 parameters. In this work we em- 
ploy the n ested restricted Monte C arlo (NRMC) method de- 
scribed in [Gregory fc Fischej|2010l to estimate the marginal 
likelihoods. Monte Carlo (MC) integration can be very in- 
efficient in exploring the whole prior parameter range be- 
cause it randomly samples the whole volume. The fraction 
of the prior volume of parameter space containing significant 
probability rapidly declines as the number of dimensions in- 
crease. For example, if the fractional volume with significant 
probability is 0.1 in one dimension then in 32 dimensions 
the fraction might be of order 10""^^. In restricted MC inte- 
gration (RMC) this is much less of a problem because the 
volume of parameter space sampled is greatly restricted to 
a region delineated by the outer borders of the marginal 
distributions of the parameters for the particular model. 

In nested RMC (NRMC) integration, multiple bound- 
aries are constructed based on credible regions ranging from 
30% to ^ 99%, as needed. We are then able to compute 
the contribution to the total integral from each nested in- 
terval and sum these contributions. For example, for the 
interval between the 30% and 60% credible regions, we gen- 
erate random parameter samples within the 60% region and 
reject any sample that falls within the 30% region. Using 
the remaining samples we can compute the contribution to 
the NRMC integral from that interval. 

The left panel of Figure [22] shows the contributions 
from the individual intervals for 5 repeats of the NRMC 
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Table 2. Marginal likelihood estimates, Bayes factors relative to model 4, and false alarm probabilities. The last two columns list the 
MAP estimate of the extra noise parameter, s, and the RMS residual. 



Model 


Periods 
(d) 


Marginal 
Likelihood 


Bayes factor 
nominal 


False Alarm 
Probability 


s 

(m s-i) 


RMS residual 

(m s^i) 


iVlQ 






2 X IQ— 






Q S 
y .o 


Ml 


(5.37) 


(4.221 ± 0.003) X 10""5 


1.4 X 10"" 


1.4 X 10^*2 


3.5 


3.6 


M2 


(5.37, 12.9) 


(1.94 ±0.01) X 10""^ 


6.5 X 10"* 


2.2 X 10-1'' 


2.4 


2.6 


Ma 


(5.37, 12.9,66.9) 
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Figure 22. Left panel shows the contribution of the individual nested intervals to the NRMC marginal likelihood for the 3 planet model. 
The right panel shows the integral of these contributions versus the parameter volume of the credible region. 



evaluation for the 3 planet model. The right panel shows 
the summation of the individual contributions versus the 
volume of the credible region. The credible region listed 
as 9995% is defined as follows. Let Xjjgg and Xlqq corre- 
spond to the upper and lower boundaries of the 99% cred- 
ible region, respectively, for any of the parameters. Simi- 
larly, Xj795 and X_L95 are the upper and lower boundaries of 
the 95% credible region for the parameter. Then Xygggs = 
Xu99 + {Xugg — Xusb) and XL9995 = Xl99 + {Xlsb — ^lqs). 
Similarly, Xj799g4 = X^/gg -I- (X^/gg — Xuha)- 

The NRMC method is expected to underestimate the 
marginal likelihood in higher dimensions and this under- 
estimate is expected to become worse the larger the num- 
ber of model pa rameters, i.e. increasing number of planets 
jCregorvl [20073 ). When we conclude, as we do, that the 
NRMC computed odds in favor of the five planet model 
compared to the four planet model is ~ 10^ we mean that 
the true odds is > 10^. Thus the NRMC method is con- 
servative. One indication of the break down of the NRMC 
method is the increased spread in the results for repeated 
evaluations. 

We can readily convert the Bayes factors to a Bayesian 
False Alarm Probability (FAP) which we define in equa- 
tion [121 For example, in the context of claiming the detec- 



tion of m planets the FAPm is the probability that there are 
actually fewer than m planets, i.e., m — 1 or less. 



FAPm ~ ^ ^ (prob.of i planets) 



(12) 



If we assume a priori (absence of the data) that all mod- 
els under consideration are equally likely, then probability 
of each model is related to the Bayes factors by 

Bi4, 



(13) 



J4 



where is the maximum number of planets in the hypoth- 
esis space under consideration, and of course B44 = 1. For 
the purpose of computing FAP™ we set N = m. Given the 
Bayes factors in Table [2] and substituting into equation 112! 



FAP5 = 



(iJo4 + -B14 + B2A + B:j4 + B44,) 



10" 



(14) 



3j4 



For the 5 planet model we obtain a low FAP ^ 10"^. The 
Bayesian false alarm probabilities for 1, 2, 3, 4, 5, and 6 
planet models are given in the fourth column of Table [5] 

Table [2] gives the NMRC Marginal likelihood estimates, 
Bayes factors and false alarm probabilities for 0, 1, 2, 3, 4, 5 
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and 6 planet models which are designated Mo, ■ • ■ , Me. The 
last two columns list the maximum a posteriori (MAP) esti- 
mate of the extra noise parameter, s, and the RMS residual. 
For each model the NRMC calculation was repeated 5 times 
and the quoted errors give the spread in the results, not the 
standard deviation. The Bayes factors that appear in the 
third column are all calculated relative to model 4. 

A summary of the 5 planet model parameters and their 
uncertainties are given in Table [S] The quoted value is 
the median of the marginal probability distribution for the 
parameter in question (except eccentricity which uses the 
mode) and the error bars identify the boundaries of the 
68.3% credible region 0. The value immediately below in 
parenthesis is the MAP estimate, the value at the maxi- 
mum of the joint posterior probability distribution. It is not 
uncommon for the MAP estimate to fall close to the bor- 
ders of the credible region. In one case, the eccentricity of 
the fifth planet, the MAP estimate falls well outside the 
68.3% credible region which is one reason why we prefer to 
quote median or mode values as well. The semi-major axis 
and M sin i values are derived from the model parame ters 
assuming a stellar mass of 0.31 ± 0.02 Mq (Delfosso e t al.l 
I2OOOI ). The quoted errors on the semi-major axis and Msini 
include the uncertainty in the stellar mass. 

Although the NRMC estimate of the Bayes factor for 
the 6 planet model is much lower than for either the 4 or 
5 planet models we can still infer the orbital parameters of 
the most probable additional planetary signal in the 6 planet 
fit. The period = 34.4 ± O.ld, the eccentricity = 0.49tE!;i7, 
and the semi-major axis and Msini are (0.140 ± 0.003 au, 
2.3ir7M®). 



6 ANALYSIS OF THE HIRES DATA AND 
COMBINATION OF HIRES AND HARPS 

In this section we present results on fits to the HIRES data 
alone and the combination of HIRES and HARPS data. 
Panel (a) of Figure [23] shows the combined HIRES (grey 
points) and HARPS (black points) data for Gl 581. Panel 
(b) shows a blow-up of a portion of the mean 4 planet model 
fit compared to the data, and panel (c) shows the residuals. 
The same reference time was used for the combi ned data set 
as fo r the HARPS only data. The HIRES data l|Vogt et al.1 
I2OIOI ) consisted of 122 velocity measurements spanning a 
range of 11 years and with quoted errors ranging from 0.53 
to 4.82 m s~^. Figure [24] shows a comparison of the velocity 
differences (HIRES-HARPS) for the nearest pairs of samples 
versus the sample time difference. A cluster of 10 points with 
sample time differences between 0.2 and 0.3s is indicated by 
the ellipse. Since these time differences are small compared 
with the shortest known orbital period of 3.15 days, they 
provide an indication of the agreement between the two sets 



^ In practice, the probability density for any parameter is repre- 
sented by a finite list of values pi representing the probability in 
discrete intervals 5X. A simple way to compute the 68.3% credi- 
ble region, in the case of a marginal with a single peak, is to sort 
the Pi values in descending order and then sum the values until 
they approximate 68.3%, keeping track of the upper and lower 
boundaries of this region as the summation proceeds. 
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Figure 23. Panel (a) shows the combined HIRES (grey points) 
and HARPS (black points) data set for Gl 581. Panel (b) shows 
a blow-up of the mean 4 planet model fit compared to the data, 
and panel (c) shows the residuals. 

of measurements. The mean velocity difference for these 10 
samples is 1.8 m s~^. 

The standard deviation of the velocity differences for 
these 10 pairs is 2.37 m s"^ For comparison the mean value 
of the quoted errors for each pair added in quadrature was 
1.94 m s"^ 

6.1 Two planet fit to HIRES data 

The results of our 2 planet Kepler periodogram analysis are 
shown in Figures 1251 1261 and [27] The upper panel of Fig- 
ure [25] shows a plot of the Logio [Prior x Likelihood] versus 
FMCMC iteration for a 2 planet fit of the HIRES data. The 
lower panel shows the values of the two unknown period pa- 
rameters versus iteration number. The two starting periods 
of 5.37 and 12. 9d are shown on the left hand side of the plot 
at a negative iteration number. The median value of the ex- 
tra noise parameter s — 2.69m s~^ compared to 2.37m s~^ 
for the HARPS 2 planet fit. Figure [26] shows the two planet 
Kepler periodogram. Only values within 5 decades of the 
maximum Logio [Prior x Likelihood] are plotted but with- 
out regard to whether the values occurred before or after 
burn-in. Two prominent periods were clearly detected: 5.37 
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Table 3. Five planet model parameter estimates from HARPS analysis. 



Parameter 


planet 1 


planet 2 


planet 3 


planet 4 


planet 5 


P(d) 


(3.14977) 


5.36871;°°°^ 
(5.36866) 


12 Q27+'006 
(12.9316) 


66.85l»;i5 
(66.747) 


3991}^ 
(387.6) 


K (m s-i) 


, ot:+0.24 
1.00_Q 22 

(1.93) 


1 c;o+0.23 
iZ.OO_g 22 

(12.39) 


o 1O+0.22 
3.iS_o.24 

(3.40) 


2 43+0.31 
^■^"^-o.ai 
(2.75) 


1 0+0.4 
1-^-0.5 
(1.62) 


e 


11+0-06 

(0.197) 


o.oi5^:«?l 

(0.022) 


11+0-05 

(0.155) 


32+0-10 

'-'■■'^-0.09 

(0.38) 


21+0-" 
'J-^l-0.21 

(0.79) 


oj (deg) 


133^75 

(140) 


(-1) 


234lg 
(234) 


004+25 
(326) 


28ltioo 
(310) 


a (au) 


0.0285t-«r6 
(0.0285) 


0406+ 0003 
(0.406) 


0730+ 0016 
(0.730) 


0.2181:005 
(0.218) 


72+0-24 
'^•'^-0.24 

(0.71) 


Msini {Me) 


(1.984) 


(15.50) 


(5.63) 


°- '-0.8 
(7.38) 


6.6^2-? 
(5.14) 


Periastron 
passage 


4182.6^0.? 
(4184) 


4182ti;^ 
(4186) 


4168.91^;^ 
(4184) 


(4134) 


3803^114 
(3828) 



(JD - 2,450,000) 




Time difference (d) 

Figure 24. A comparison of the velocity differences (HIRES- 
HARPS) for the nearest pairs of samples versus the sample time 
difference. 



and 12. 9d. The second period parameter exhibited two other 
peaks but these were significantly less probable. The most 
probable of these has a period ~ 9d. 

Figure [27] shows a plot of eccentricity versus period for 
a sample of the FMCMC parameter samples for the 2 planet 
model. The dominant 5.37 and 12. 9d peaks and the weaker 
9d peak allow for low eccentricity orbits. The peak around 
300d has a high value of eccentricity typical of noise. 
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Iterations (x 10^) 
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Iterations (x 10^) 

Figure 25. The upper panel is a plot of the Logio [Prior X Likeli- 
hood] versus iteration for a sample of the FMCMC two planet fit 
of the HIRES data. The lower shows values of the two unknown 
period parameters versus iteration number. The two starting pe- 
riods of 5.37 and 12. 9d are shown on the left hand side of the plot 
at a negative iteration number. 



6.2 Three planet fit to HIRES data 

Two three planet runs were carried out on the HIRES 
only data starting with the best periods (5.37,12.9, 66. 9d) 
found from the HARPS analysis but only the 5.37d period 
(largest amplitude) was successfully detected. The best of 
these two runs detected three dominant periods of 5.37, 8.99 
and~ 300d. A much weaker peak was found at a period of 



12. 9d. The HIRES fit extra noise parameter was s = 2.2 m 
s~^ compared to the HARPS 3 planet fit where s — 1.7 m 
s~^. Figure [29] shows a plot of the three planet Kepler pe- 
riodogram. Only values within 5 decades of the maximum 
Logio [Prior x Likelihood] are plotted but without regard to 
whether the values occurred before or after burn-in. Three 
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Figure 26. A plot of the 2 period parameter values versus a nor- 
malized value of Logio [Prior X Likelihood] , the 2 planet Kepler 
periodogram of the HIRES data. 
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Figure 27. A plot of eccentricity versus period for the 2 planet 
FMCMC fit of the HIRES data. 



prominent periods were clearly detected: 5.37, 9, and 300d. 
The second period parameter (shown in grey) exhibited a 
second much less probable peak at 12. 9d and the third pe- 
riod parameter (black) exhibited many weak peaks. 

Figure [30] shows a plot of eccentricity versus period for 
a sample of the FMCMC parameter samples for the 3 planet 
model. The 5.37 and 9d peaks exhibit low eccentricity orbits. 
The peak around 300d has a high value of eccentricity. 

Figure [31] shows a blow-up of the above eccentricity 
versus period plot in the vicinity of the 300d peak which is 
dominated by two high eccentricity features typical of noise. 
We conclude that there is no clear evidence for a third period 
in the HIRES data alone and suspect that the presence of 
the strong high eccentricity 300d complex may contribute 
to the dominance of the 9d period over the 12.9 day period 
found in the 2 planet fit. 

For both the 2 planet and 3 planet fits the extra noise 
parameter is larger for the HIRES data than the HARPS 
data. One possibility is that the quoted HIRES errors have 
been systematically under-estimated which we can model 
by an extra Gaussian noise term added in quadrature to 
the quoted HIRES errors with a. a — dsHiRES- We can 
obtain a crude estimate of dsniRES from the HIRES and 
HARPS s parameter values for the 2 planet case where 
both analyses yielded the same 2 periods. The result is 
dsHiRES = V2.692 - 2.372 = 1.3 m s"\ This suggest that 
in the analysis of the combined data set we should include 
an extra Gaussian noise term added in quadrature to the 
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Figure 28. The upper panel is a plot of the Logio [Prior X Like- 
lihood] versus iteration for a sample of the FMCMC three planet 
fit of the HIRES data. The lower shows the values of the three 
unknown period parameters versus iteration number. The start- 
ing periods of 5.37, 12.9 and 66. 9d are shown on the left hand 
side of the plot at a negative iteration number. 




5.37 9 12.< 



300 1000 

Period (d) 



Figure 29. A plot of the 3 period parameter values versus a nor- 
malized value of Logio [Prior x Likelihood], the 3 planet Kepler 
periodogram of the HIRES data. 



quoted HIRES errors with the a, labeled dsniRES, as an ad- 
ditional unknown parameter. 



6.3 Two planet fit to the combined 
HIRES/HARPS data 

Based on the above results we decided to use the follow- 
ing noise model for the j^^ data point for the combined 
HIRES/HARPS two planet analysis. 



.2 

' jquotcd 



- ^ ' '^iquotcd + '^•^IRES 



(15) 



(16) 



18 P. C. Gregory 



\ 




■ 

i 

. i 

\ 





3. 5.37 9 12.9 30 100 300 1000 3000 10000 30000 



Period (d) 

Figure 30. A plot of eccentricity versus period for the 3 planet 
FMCMC fit of the HIRES data. 
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Figure 31. A blow-up of the above eccentricity versus period 
plot in the vicinity of the 300d peak. 



The best two planet orbital parameters were employed as 
start coordinates for the combined HIRES/HARPS analy- 
sis. We also incorporated an additional unknown parameter 
dc to allow for a possible difference in the constant velocity 
offsets of the HIRES and HARPS data. Based on our anal- 
ysis of Figure [24l our best estimate of dc ~ 1.8m s~^. We 
assumed a Gaussian prior with zero mean and a = 3m s^^. 

The two planet FMCMC fit of the HIRES /HARPS data 
confirmed the 5.37 and 12. 9d periods. Figure [32] shows a plot 
of a subset of the FMCMC parameter marginal distributions 
for the 2 planet fit of the data after filtering out the post 
burn-in FMCMC iterations that correspond to the 2 domi- 
nant period peaks at 5.37 and 12. 9d. The bottom row shows 
the marginals for dc, dsHARPS, and s. The maximum a pos- 
terior and median values found for dc are 1.65 and 1.67m 
, respectively, i.e., very close to the crude estimate of 
1.8m s~^ made earlier. 

The upper panel of Figure [33] shows the marginal distri- 
bution for the unknown dsniRES parameter in the 2 planet 
FMCMC fit of the combined HIRES/HARPS data. The 
black curve in the lower panel shows the marginal dis- 
tribution for the common extra noise parameter s in the 
HIRES/HARPS fit. The light grey curve is the same quan- 
tity for the two planet HARPS only fit. Clearly there is very 
good agreement between the two s parameter estimates. 
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Figure 32. A plot of a subset of the FMCMC parameter marginal 
distributions for a 2 planet fit of the combined HIRES/HARPS 
data. 
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Figure 33. The upper panel shows the marginal distribution for 
the unknown dsniRES parameter in the 2 planet FMCMC fit of 
the combined HIRES/HARPS data. The black curve in the lower 
panel shows the marginal distribution for the common extra noise 
parameter s in the HIRES/HARPS fit. The light grey curve is the 
same quantity for the two planet HARPS only fit. 



6.4 Three planet fit to the combined 
HIRES/HARPS data 

In the above two planet fit to the combined HIRES/HARPS 
we found that the marginal distribution for the extra noise 
parameter s agreed closely with that obtained from the two 
planet HARPS only data. For the three planet combined 
analysis we decided to fix the s parameter to a value of 
1.9m s~^, i.e., the value obtained from the HARPS alone 3 
planet fit. We decided to use the following noise model for 
the J*'' data point for the three planet analysis. 

'SHARPS, = y^<,ot.d, + 1-9^ (17) 
'^HIRESj = l/o-quotodj + rf4lRES + l-S'' (18) 



Bayesian Re-analysis of the Gliese 581 Exoplanet System 19 




Figure 34. A plot of the 3 period parameter values versus a 
normalized value of Logio [Prior X Likelihood] for the 3 planet 
FMCMC Kepler fit of the combined HIRES/HARPS data. 
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Figure 35. A plot of eccentricity versus period for the 3 planet 
FMCMC fit of the combined HIRES/HARPS data. 



Again the best three planet orbital parameters from the 
HARPS analysis were employed as start coordinates for the 
combined HIRES/HARPS analysis. All 3 period parameters 
were free to roam within a search range extending from l.ld 
to 10 X the data duration. 

Figure l34l shows plot of a sample of the FMCMC three 
planet periodogram. Only values within 5 decades of the 
maximum Logio [Prior x Likelihood] are plotted but with- 
out regard to whether the values occurred before or af- 
ter burn-in. Three prominent periods were clearly detected: 
5.37, 12.9, and 66. 9d. The third period parameter exhibited 
other peaks but these were significantly less probable, the 
one at 82d is consistent with being a one year alias of the 
66. 9d period. 

Figure [SS] shows a plot of eccentricity versus period for 
a sample of the FMCMC parameter samples for the 3 planet 
model. The dominant 5.37, 12.9, and 66. 9d peaks allow for 
low eccentricity orbits. 

Figure [35] shows a plot of a subset of the FMCMC pa- 
rameter marginal distributions for the 3 planet fit of the 
data after filtering out the post burn-in FMCMC iterations 
that correspond to the 3 dominant period peaks at 5.37, 
12.9, and 66. 9d. The bottom row shows the marginals for 
dc, dsHiRES. The dsniRBS is more accurately defined than 
in the two planet analysis. The maximum a posterior and 
median values found for dsniRBS are 1.61 and 1.69m s~^, 
respectively. 




dc (m s ^ dSHiHEs (ms\ 

Figure 36. A plot of a subset of the FMCMC parameter marginal 
distributions for a 3 planet fit of the combined HIRES/HARPS 
data. 



6.5 Four planet fit to the combined 
HIRES/HARPS data 

For the four planet fit we reverted to the noise model used 
for the two planet analysis. 



CHARPSi 



CHIRESi ~ 



quotcdj 



+ ds? 



(19) 



(20) 



The best four planet orbital parameters from the HARPS 
only analysis were employed as start coordinates for the 
combined HIRES/HARPS analysis. As before we incorpo- 
rated the additional unknown parameter dc to allow for 
a possible difference in the constant velocity offsets of the 
HIRES and HARPS data. 

The four planet Kepler periodogram found the four 
starting periods of 3.15, 5.37, 12.9, and 66. 9d and no other 
peaks. Figure [571 shows the marginal posterior densities of a 
subset of the parameters. 

Table [4] shows a comparison of the estimates of the dc, 
dsHiRES, and s parameters from the 2, 3, and 4 planet fits 
to the combined HIRES/HARPS data set. The parameter 
value listed is the median of the marginal probability dis- 
tribution for the parameter in question and the error bars 
identify the boundaries of the 68.3% credible region. The 
value immediately below in parenthesis is the MAP esti- 
mate, the value at the maximum of the joint posterior prob- 
ability distribution. As expected the extra noise parameter 
s decreases with the number of planets fit. The estimates 
of the offset parameter, dc, and the HIRES extra noise pa- 
rameter, dsHiRES, agree within the quoted uncertainties and 
these uncertainties are smallest for the 4 planet model fit. 

We carried out an analysis of the 4 planet normalized 
fit residuals (fij)- We define ri. in equ. 1211 



ri. 



rii ^ 



residuali 



^ ^crror2 + A^\i^^^.^ sf ' 



(21) 



where j is an index for the combined HIRES/HARPS data 
set. rii is the number of post burn-in FMCMC equilib- 
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Figure 38. A histogram of tlie values compared to a suitably 
normalized Gaussian with zero mean and standard deviation = 1. 



Figure 37. A plot of a subset of the FMCMC parameter marginal 
distributions for a 4 planet fit of the combined HIRES/HARPS 
data. 



Table 4. Bayesian estimates of parameters dc, rfsniRES > and s 
from the combined HIRES/HARPS data set for the 2, 3, and 4 
planet fits. 



Parameter 


2 planet 


3 planet 


4 planet 


dc 

(m s-i) 


^■"'-0.35 
(1.64) 


1 c:c + 0.25 
J^-00_0.44 
(1.65) 


1 45+0-32 
(1.53) 


(iSHIRES 
(m s~i) 


(0.75) 


1 fiq+O^'S 
J^-O»_0.32 

(1.61) 


1 04+0.35 

(1.85) 


s 

(m s~-^) 


o qQ+0.25 
^••^^-0.11 

(2.32) 


1.9 

fixed 


1.34+;}^ 
(1.45) 
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Figure 40. A plot of the period parameter versus a normalized 
value of Logio [Prior X Likelihood] for the 1 planet FMCMC Ke- 
pler fit to the HIRES residuals after removing the 3.15, 5.37, and 
12.9 day orbits.. 



rium samples used in computing the mean value of resid- 
ual/ (effective noise a) for each measurement (typically we 
use rii ~ 200). The effective noise a is given by 

effective noise a — ^ error? -|- ds^j^^g. -|- s?, (22) 

which is the quadrature sum of the quoted error, dsniRESi > 
and the extra noise term Si. Of course, dsniRESi only con- 
tributes to the effective noise cr of the HIRES data values. 

Figure l38] shows a histogram of the n^. values compared 
to a suitably normalized Gaussian with zero mean and stan- 
dard deviation = 1. Within the uncertainties ri . is consis- 
tent with a Gaussian distribution. The reduced x of the ri- 
is given by equ. [23l 

X? = V rl = 1.008, (23) 

nt — Up ' 
j=i 

where nt — the total number of combined HIRES/HARPS 
data points and Up = 24 is the number of fit parameters. 

The four period phase plots are shown in Figure 1391 
The top left panel shows the data and model fit versus 



3.15d orbital phase after removing the effects of the three 
other orbital periods. The upper and lower solid curves are 
the mean FMCMC model fit ±1 standard deviation. The 
HIRES data points are shown in grey and the error bars are 
the quoted errors added in quadrature with our median esti- 
mate of dsHiRESi ~ 1.84m s~^. The other panels correspond 
to phase plot for the other three periods. Examination of the 
66. 9d period phase plot indicates that the HIRES data do 
not lend much support to the 66.9 day period. In fact around 
a phase of 0.8 there are a series of ~ 11 points that exhibit 
a systematic trend away from the mean light curve. To ex- 
amine this further we carried out a one planet FMCMC fit 
to the HIRES residuals after subtracting off the 3.15, 5.37, 
and 12. 9d mean orbits. Figure HQ] shows the one planet Ke- 
pler periodogram of these residuals. Three prominent peaks 
are present at 26.3, 65.6, and 73d along with many minor 
peaks. The thin solid vertical lines highlight these peaks. 
The strongest peak has a period of 73d which explains the 
absence of good support for a 66. 9d period in Figure l39l 
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Figure 39. Phase plots for the 4 most probable periods derived from the combined HIRES/HARPS data. The top left panel shows the 
data and model fit versus 3.15 day orbital phase after removing the effects of the five other orbital periods. The upper and lower curves 
are the mean FMCMC model fit ±1 standard deviation. The other three panels correspond to phase plot for the other three periods. 
The HARPS data points are black and the HIRES grey. 
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Figure 41. A plot of the 5 period parameters versus a normal- 
ized value of Logio [Prior X Likelihood] for the 5 planet FMCMC 
Kepler fit to the combined HIRE/SHARPS data. 



6.6 Five planet fit to the combined 
HIRES/HARPS data 

Figures HT] and show tlie results of a five planet fit to the 
combined HIRES/HARPS data set. The evidence for a ~ 
400d period is more confused than the IfARPS only results 
shown in Figures [H [13 The combined HIRES/HARPS re- 
sults show a dominant high eccentricity peak with a period 
of 472d. Again, weak high eccentricity peaks are a common 
characteristic of noise. 



T DISCUSSION 

Our analysis of the combined HIRES/HARPS data set ar- 
gues strongly that the HIRES errors have been systemati- 
cally underestimated. If we model this by the additive extra 
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Figure 42. A plot of eccentricity versus period for the 5 planet 
FMCMC fit to the combined HIRE/SHARPS data. 



noise term ^shires , as in equ. 1241 we conclude dsHiRES 



1.84 



+0.35 



An alternative possibility that we can test is that 



the HIRES data are systematically too low by a common 
factor. To check this out we re-did the four planet analysis 
of the combined HIRES/HARPS data using the following 
noise model. 



f7"HIRESi 



\J (&HIRES X (Tquotcdj)^ 



+ S2 



(24) 



We assumed a uniform prior for the unknown parameter 
&HIRES in the range 0.5 to 4.0. The results were qualita- 
tively very similar to the results obtained with dsniRES noise 
term. Following similar calculations to those outlined in Sec- 
tion 16.51 we computed the reduced of the residuals di- 
vided by the total effective noise cr and obtained a value of 
1.009, indistinguishable from the reduced obtained us- 
ing the dsHiRES parameter. Figure [43] shows the marginal 
probability distributions for 6hires and s parameters. The 
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Figure 43. A plot of the marginal probability distributions for 
the bniRES s parameters for a 4 planet FMCMC fit to the 
combined HIRE/SHARPS data using the noise model of equ. 1241 
The dashed curve in the lower panel is the marginal for s from 
the 4 planet fit to the HARPS only data. 



dashed curve in the lower panel is the marginal for s from 
the 4 planet fit to the HARPS only data. Employing the 
Shires parameterization results in a significantly larger es- 
timate in the s term than was required by the HARPS only 
analysis or the dsHiRES parameterization of the combined 
HIRES/HARPS data. 

Although the two different parameterizations of the 
HIRES extra noise lead to similar values of the reduced 
we can gain additional insight about which parameterization 
is better from a more microscopic exploration by binning the 
4 planet fit residuals/ (effective noise) and examining the 
values of the individual bins. If we add two random vari- 
ables one of which has v\ degrees of freedom and the other 
has V2 degrees of freedom then thei r sum will be with 
v\ -\- V2 degrees of freedom, e.g., see lGregorvll2005al p. 144. 
Similarly, we can take the quantity ri- of equ. [23] which is 
with {nt — n-p) degrees of freedom and divide the ri- values 
into k bins. Suppose the k^^ bin has ri - values then the 
of these values will be given by equ. 1251 



2 

Xk 



^x(nt - Up) ^ 



(25) 



which is X with — x (nt 



degrees of freedom. 



Figure 33] shows a plot of the x^ statistic versus the 
original quoted errors that have been binned into 0.5m s~ 
bins for the 4 planet fit to the combined HIRES/HARPS 
data using the dsniRES parameterization. Figure [45] shows a 
similar plot for the 6hires parameterization. Ideally the x^ 
distribution would be flat with a value of 1.0. It is clear that 
the dsHiRES parameterization achieves a flatter distribution 
than the &hires parameterization. Further the noise model 
involving dsniRES leads to a jitter noise estimate s which is 
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Figure 44. A plot of the x| statistic from equ. I25l versus binned 
values of the quoted errors for the 4 planet fit of the combined 
HIRES/HARPS data using the rfsniRES parameterization. The 
solid and dashed histograms show the binned HIRES data and 
binned HARPS data, respectively. 
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Figure 45. A plot of the statistic from equ. I25l versus binned 
values of the quoted errors for the 4 planet fit of the combined 
HIRES/HARPS data using the 6hires parameterization. The 
solid and dashed histograms show the binned HIRES data and 
binned HARPS data, respectively. 



much closer to the value expected from the HARPS only 4 
planet analysis. 

Our HARPS only analysis provides evidence for 5 plan- 
etary signals while incorporating the HIRES data seems to 
degrade the evidence for both a 66.9 and ~ 400d periods 
even when we allow for an extra HIRES noise term of order 
1.8m s"'^. We suspect this extra noise term may arise from 
an as yet unidentified systematics which may be the reason 
for the degradation in the quality of fits beyond three plan- 
ets. Alternatively, could unidentified HARPS systematics be 
responsible for the extra periods evident in their lower noise 
data? 



8 CONCLUSIONS 

A Bayesian re-analysis of published HARPS and HIRES 
precision radial velocity data for Gl 581 was carried out 
with a multi-planet Kepler periodogram (from 1 to 6 plan- 
ets) based on our fusion Markov chain Monte Carlo algo- 
rithm. In all cases the analysis included an unknown pa- 
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rameterized stellar jitter noise term. For the HARPS data 
set the most probable number of planetary signals detected 
is 5. The Bayesian false alarm probability for the 5 planet 
model is 0.01. These include the 3.15, 5.37, 12.9, 66. 9d pe- 
riods reported previously plus a 399ti6d period. The or- 
bits of 4 of the 5 planets are consistent with low eccen- 
tricity orbits, the exception being the 66. 9d orbit where 
e = 0.3310 50. semi-major axis and Msini of the 5 plan- 
ets are (0.0285 ± 0.0006 an, I.9I+O25M0), (0.0406 ± 0.0009 
an, 15.7t^jM0), (0.0730 ± 0.0016 au, 5.29 ± O.43M0), 
(0.218 ± 0.005 au, 6.7 ± O.8M0), and (0.72 ± 0.24 au, 
6.6l2^°M0), respe ctively. 

In light of the lVogt et al.l (|2010l ) report of a sixth com- 
panion with a period of 36. 6d, we carried out a 6 planet fit 
to the HARPS data which detected multiple period possibil- 
ities. The strongest of these, with a period — 34.4±0.1d and 
eccentricity of 0.49j;Q j7, had a peak Logio [Prior x Likeli- 
hood] 100 times larger than the others. The inferred semi- 
major axis and Msini are (0.140 ± 0.003 au, 2.3+°-^M(^). 
The Bayesian false alarm probability for the six planet model 
is extremely large 0.999978 so we are unable to support any 
claim for a sixth companion on the basis of the current data. 

The analysis of the HIRES data set yielded a reliable 
detection of only the strongest 5.37 and 12.9 day periods. 
The analysis of the combined HIRES/HARPS data again 
only reliably detected the 5.37 and 12. 9d periods. Detection 
of 4 planetary signals with periods of 3.15, 5.37, 12.9, and 
66. 9d was only achieved by including an additional unknown 
but parameterized Gaussian error term added in quadrature 
to the HIRES quoted errors. The marginal probability den- 
sity of the sigma for this additional HIRES Gaussian noise 
term has a well defined peak at 1.84tQ'33m s~^. Phase plots 
indicate that incorporating the HIRES data seems to de- 
grade the evidence for the 66.9 and ~ 400d periods even 
when we allow for the extra HIRES noise term. We suspect 
this extra noise term may arise from unidentified systematics 
which may be the reason for the degradation in the quality 
of fits beyond three planets. Alternatively, could unidenti- 
fied HARPS systematics be responsible for the extra periods 
evident in their lower noise data? Independent experimental 
confirmation of the HARPS low noise level results would be 
very desirable. 
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